Intervalley Coupling for Interface-Bound Electrons in Silicon: 
An Effective Mass Study 
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Orbital degeneracy of the electronic conduction band edge in silicon is a potential roadblock to the 
storage and manipulation of quantum information involving the electronic spin degree of freedom 
in this host material. This difficulty may be mitigated near an interface between Si and a barrier 
material, where intervalley scattering may couple states in the conduction ground state, leading to 
non-degenerate orbital ground and first excited states. The level splitting is experimentally found 
to have a strong sample dependence, varying by orders of magnitude for different interfaces and 
samples. The basic physical mechanisms leading to such coupling in different systems are addressed 
here. We expand our recent study based on an effective mass approach, incorporating the full plane 
wave (PW) expansions of the Bloch functions at the conduction band minima. Physical insights 
emerge naturally from a simple Si/barrier model. In particular, we present a clear comparison 
between ours and different approximations and formalisms adopted in the literature, and establish 
the applicability of these approximations in different physical scenarios. 
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I. INTRODUCTION 

Electronic spins in Si are promising candidates for 
qubits due to their naturally long coherence times»i An 
important challenge for using electronic spins as qubits 
in silicon devices is to assure that the low energy physics 
is ruled solely by the 2-level spin degree of freedom^r— 
In bulk Si crystal, the conduction band lower edge is 
six-fold degenerate. Any superposition of the six Bloch 
states associated with the minima in k along the ±x, ±y 
and ±z crystallographic directions is also an eigenstate 
of the crystalline hamiltonian, so that the orbital state of 
a free conduction electron at the band edge is normally 
not defined. 

Several quantum computer architectures under inves- 
tigation involve manipulation of the electronic spin at 
an interface between Si and some barrier material, most 
commonly SiGe alloy o^^d^ and SiO^r^^ barriers. The 
confining electric field (generated by external electro- 
static gates and/or depletion layers) generates a quasi- 
triangular potential well at the interface. Assuming the 
interface to be perpendicular to the z direction [that is, 
a (001) interface], the ground state energy of such tri- 
angular potential well depends on the effective mass in 
the z direction. The effective masses at the conduction 
band minima of Si arc anisotropic, with the longitudi- 
nal effective mass more than 4 times larger than the 
transversal. This shifts the minima along the x and 
y directions well above the minima in the z direction, 
breaking the six-fold degeneracy into a 2-fold degenerate 
ground state and a 4-fold degenerate excited statcj ^^d'' 
The splitting is further enhanced if tensile strain is ap- 



plied to the Si crystal {e.g. in Si quantum wells grown 
over relaxed SiGe substrates) The degeneracy in the 
two-dimensional {fc^, k-z = —k^} subspace is lifted in 
the presence of a sufficiently singular perturbation poten- 
tial, such as a Si/barrier interface.— Experimental values 
of the valley splitting in interfaces have been reported in 
the 0.1 to 1 meV range.— A peculiar result was re- 
ported by Takashina et al.^ who measured a "giant" 
splitting of 23 meV on a Si/Si02 interface in a SIMOX 
(separation by implantation of oxygen) structure. 

In the presence of orbital degeneracy, electron manip- 
ulations relying on the Pauli's exclusion principle, such 
as Heisenberg exchange coupling^^ and spin blockade^^, 
may become unreliable since the qubits Hilbert Space is 
now spanned by the valley as well as the spin degrees 
of freedom. For an electron spin qubit confined in a 
Si quantum dot, reliable knowledge of the interface in- 
duced valley splitting is crucialJ^ For example, for a sin- 
gle electron spin qubit if the valley splitting is smaller 
than the reservoir thermal energy, or valley splittings are 
different across two quantum dots, exchange gates can- 
not be performed as was originally designed since two- 
valley two-electron singlet and triplet states (where the 
two electrons are in different valleys) are not exchange- 
split When two-electron singlet and triplet states are 
used to encode a logical qubit, reliable initialization be- 
comes impossible if valley splitting in a quantum dot is 
unknown or known to be small (compared to reservoir 
temperature energy scale) »21 In short, clear knowledge of 
a large valley splitting in a Si quantum dot is imperative 
to assure the feasibility of an electron spin qubit. 

Early theoretical studies on Si valley splitting in the 
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framework of the EMA were performed more than 30 
years agOi ^^i^^ The relevance of the periodic part of the 
Si bulk Bloch states (leading to the so-called Umklapp 
processes) to this scattering became clear in the 70'Sf21 
However, its inclusion combined with the EMA formal- 
ism leads to a puzzling dependence of calculated physical 
properties on the particular position of the barrier within 
the range of separation between atomic layers, an arti- 
fact already found in previous works^^ which is discussed 
in detail and clarified in Ref. [l^. These early theoreti- 
cal studies do not discuss the strong sample dependence 
observed experimentally. 

Valley mixing has been extensively studied in the liter- 
ature^Si"^ regarding GaAs/AlAs or GaAs/GaAlAs inter- 
faces and superlattices. The most relevant valley mixing 
in this context involves the X-valleys of the electron con- 
duction band in the AlAs layer, which is the global band 
minimum in AlAs and a local minimum in GaAs, and the 
F-valley in the GaAs layer. This mixing involves sepa- 
rated spatial layers. Higher F and X states in AlAs and 
GaAs also intervene in the mixing^^. In Si, mixing oc- 
curs between the crystalline momenta kz and —k^ in the 
same spatial layer, that is, in the Si slab, and is due to 
the barrier potential alone. The theory and experimen- 
tal consequences of GaAs/ AlAs valley mixing are thus 
different and not transferrable to Si. 

Performing a fully ah initio treatment is not realistic 
in the study of valley coupling in Si due to limitations 
both on the length and on the energy scales. Interval- 
ley splittings of the order of tenths of meV would not 
be accurately resolved within the density functional the- 
ory (DFT) approach based on current computational re- 
sources. Also, the electronic states under study spread 
over several lattice parameters, and simulation of large 
superccUs with appropriate description of the band gap 
(through Hcdin's GW schemci^ for example), again in- 
volve numerical computations beyond currently available 
capability. 

More recent investigations often employ approaches 
with atomistic ingredients, whether based completely on 
the tight-binding (TB) method or on the hybrid of TB 
and EMAi^"— These atomistic methods allow treat- 
ment of disorder effects directly. For example, TB and 
EMA-I-TB calculations conclude that for tilted Si/SiGe 
quantum wells, alloy disorder and interfacial step disor- 
der must be included to obtain finite valley splittings 
Such sample-dependent results arc consistent with the 
observed variation of this coupling as measured in differ- 
ent interfaces. On the other hand, the detailed descrip- 
tion of disorder comes at the expense of generality and 
prevents analytical insights since these methods require 
numerical treatment (with exceptions, some of which are 
discussed in Sec. IIVB]) . 

We have recently performedi^ a study of the valley 
splitting problem that leads to the identification of rele- 
vant physical mechanisms underlying the intervalley cou- 
pling. Our EMA model incorporates the atomistic Bloch 
functions obtained from ab initio calculations in connec- 



tion with an envelope function obtained from the single 
valley EMA equation. That study focused on the rele- 
vance of the barrier height and the interface width. 

In the present work, we discuss the ab initio results 
which are required in implementing our approach. In par- 
ticular, wc provide a detailed roadmap to the connection 
between effective mass and and the ab initio wavefunc- 
tion, deriving the formalism and discussing the physics 
and applicability of our methodology. Further results are 
presented and discussed, including a more complete anal- 
ysis of the role of the confining electric field perpendicular 
to the interface plane. The plane wave (PW) expansion 
of the periodic part of each Bloch function is explicitly 
given, and contributions from reciprocal lattice vectors 
in the PW expansion are analyzed, identifying the most 
relevant Umklapp processes. We also discuss the present 
model in the context of existing studies and compare our 
results with previous ones whenever warranted. The in- 
terface is modeled here within EMA by a finite height 
step potential. Other interface profiles were already con- 
sidered in our previous studyJ^ 

In section |TT] we briefly review the model, emphasizing 
the assumptions involved in the derivation of the formal- 
ism. Wc proceed to numerically calculate the coupling 
under various electric fields and conduction band offsets 
in section [mi These results allow comparison of the val- 
ley coupling in different experimental conditions, and to 
discuss possible mechanisms contributing to the giant 
splitting reported in Ref. [21I . Comparison of our theo- 
retical approach with others in the literature is given in 
section IIV[ where we also evaluate some of the approxi- 
mations adopted in previous works and shed light in some 
theoretical questions, such as the range of applicability 
of various EMA treatments. Finally, our conclusions are 
presented in section IVl 

II. THEORETICAL BACKGROUND 

We consider a single electron at the bottom of the Si 
conduction band, near a (001) Si/barrier interface. As- 
suming translational symmetry in the xy plane, we model 
the barrier material as an effective potential mimicking 
the conduction band offset along the z direction. This 
model addresses only the position in energy of the bot- 
tom of the conduction band, disregarding the detailed 
electronic structure in the transition region between Si 
and the barrier material, as we discuss below. 

The Hamiltonian for the conduction band electron is 
assumed to be of the form 



H = Ho + Uiz)-e—z, (1) 
e(z) 

where Hq is the unperturbed bulk Si Hamiltonian and 
U{z) is the barrier potential. An electric field F/e^z) 
along the z direction keeps the electron close to the in- 
terface. The dielectric function e(z) changes from the 
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bulk Si value to the bulk barrier material value and, in a 
more accurate analysis, could include many-body effects. 

A sequence of approximations are involved in obtain- 
ing the EMA equation for Si heterostructures. The first 
assumption involves the expansion of the electronic wave- 
function in terms of the periodic part of Bloch states 
around the conduction band minima, assumed to be the 
same for both the well semiconductor and the barrier 
matcrial4^ The envelope function approach is better jus- 
tified if the wavevectors at the minima of the two ma- 
terials are relatively near relative to the Brillouin zone 
dimensions. 

We further assume their effective masses to be the 
same, although it is possible within EMA to account for 
different effective masses in heterostructures.— We avoid 
specifying the barrier material and systematically inves- 
tigate the effects of the barrier height alone. 

Finally, the crystalline structure is assumed to be pre- 
served across the interface. In a MOSFET geometry, 
Si02 grown over Si is typically amorphous, as differ- 
ent metastable crystalline phases coexist, with crystallo- 
graphic directions that do not necessarily match the di- 
rections of the Si slab4^ Moreover, the most energetically 
favorable crystalline phases of Si02 have a single conduc- 
tion band minimum at the F point4^ On the other hand, 
Si(i_2.)Ge2: alloys at low concentrations of Ge {x < 0.3) 
present conduction band minima at the same position in 
the Brillouin zone as pure Si crystal;^ and the effective 
mass remains almost unchanged for the Gc concentration 
in the alloy up to x ^ 0.3. For samples under experimen- 
tal investigation)^ Si is epitaxially grown over a relaxed 
SiGe substrate, so that the crystallographic directions 
match. Although the conduction band states are not ex- 
actly the same for these two materials, they are expected 
to be very similar. Thus, for Si/SiGe heterostructures 
the EMA assumptions are better justified than for MOS- 
FETs. 

However, given the large conduction band offset be- 
tween Si and Si02 a very small penetration of the enve- 
lope function into the barrier material is expected. De- 
tailed simulation of the electronic structure of Si/Si02 
interfaces is possible, but does not lead to general re- 
sults, since different nanofabrication methods and small 
variations of the growth parameters lead to very differ- 
ent interface morphologies We adopt the EMA ap- 
proach in the case of Si/Si02 interface bearing in mind 
that this approximation could lead to quantitative inac- 
curacies that should be estimated using some other ap- 
proach. 

Within the above assumptions, one obtains the enve- 
lope function ^'(z) from the single-valley effective mass 
equation^ 



2toz dz"^ 



U{z)- 



F 
e-r-rz 

e{z) 



■^{z)^E^{z), (2) 



where is the longitudinal effective mass for Si. The 
electronic eigenstates of Eq. ([1]) bound to the interface 
are obtained from the sin gle- valley EMA^ as 0^(r) = 



\E'(z)e ^'^M^(r) where = ifegz are the Bloch wave 
vectors of the conduction band minima (fep ~ 0.84 x 
27r/ao). It is convenient to perform a PW expansion of 
the periodic functions u±{y)^ leading to the complete 
wavefunctions 



(/)±(r) = ^{z)e 



±ikoz \ ^ 



G 



(3) 



where {G} are reciprocal lattice vectors. The PW expan- 
sion in equation ([3]) , originally explored in connection to 
Ref. m, is useful in several other contexts. It was ob- 
tained from ab initio Density Functional Theory (DFT) 
calculations^ performed with the ABINIT code4^ From 
a DFT perspective, the electronic correlations for the 
bulk Bloch states are described in the Local Density Ap- 
proximation (LDA)i^i^ The exchange-correlation po- 
tential parameterized by Perdew and Zunger— from 
Ceperley- Alder quantum Monte-Carlo results for the ho- 
mogeneous electron gaa^ was adopted. The interactions 
between valence electrons and ions are described by the 
ab-initio, norm-conserving pseudopotcntials of TrouUier- 
Martins^Sst generated with the FIII98PP code^ These 
approximations significantly speed up the computation 
of the conduction band structure, with a PW expansion 
of the wavefunctions including terms up to 16 Ry, that is, 
corresponding to 290 plane waves for each fc^, with virtu- 
ally no computational effort. The calculated equilibrium 
lattice constant of Si at oq ~ 5.41 A and the conduction- 
band minima at ko — 0.844(27r/ao) are in close agreement 
with experimental resultsi^ These values are used in the 
calculations presented below. 

We find that over 90% of the spectral weight of the 
PW expansion in equation ([3]) comes from the five points 
in the BCC reciprocal lattice which are nearest to each 
conduction band minimum fc^. The valley coupling 



Vvo = {^+\H\<I>^ 



(4) 



is the key quantity leading to the valley splitting A = 
2|VVo|,i^ Therefore, a preliminary estimate of the ±z 
intervalley coupling could involve 9 PWs (the F point is 
a common nearest neighbor for both at the minima). 
The coefficients such that |c^(G)p > 10~^, accounting 
for 99.5 % of the spectral weight, are explicitly given for 
/i = z in Table [U coefficients for all band minima may be 
obtained from those by symmetry (see Table caption). 
From the above expansion, eq. ^ reads 



Vvo^ E c*+{G)c.{G')S{G,-G'J6{Gy-G'y)I{G,,G',), 

(5) 



G,G' 



where the last term stands for the integral 



/-too 
|vI/(z)|2e'Q- 
-oo 

with Q = Gz - - 2fco- 



U{z)~^z 



dz (6) 
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G 


Re[c+,(G)] 


Im[c+.(G)] 




(1-1-1) 


-0.3131 


-0.3131 


0.1961 


(-1 1 -1) 


-0.3131 


-0.3131 


0.1961 


(11-1) 


-0.3131 


0.3131 


0.1960 


(-1 -1 -1) 


-0.3131 


0.3131 


0.1960 


( 0) 


0.3428 


-0.0000 


0.1175 


(2 0-2) 


-0.0986 


0.0000 


0.0097 


(0 2-2) 


-0.0986 


0.0000 


0.0097 


(-2 -2) 


-0.0986 


0.0000 


0.0097 


( -2 -2) 


-0.0986 


0.0000 


0.0097 


(1-11) 


0.0695 


-0.0695 


0.0097 


(-1 1 1) 


0.0695 


-0.0695 


0.0097 


(111) 


0.0695 


0.0695 


0.0097 


(-1 -1 1) 


0.0695 


0.0695 


0.0097 


(-2 2 -2) 


-0.0000 


-0.0451 


0.0020 


( 2 -2 -2) 


-0.0000 


-0.0451 


0.0020 


(-2 -2 -2) 


0.0000 


0.0451 


0.0020 


(2 2-2) 


0.0000 


0.0451 


0.0020 


( 2 2) 


0.0387 


-0.0000 


0.0015 


( 2 2) 


0.0387 


-0.0000 


0.0015 


(0-2 2) 


0.0387 


-0.0000 


0.0015 


(-2 2) 


0.0387 


-0.0000 


0.0015 


(0 0-4) 


0.0186 


-0.0000 


0.0003 


(-1 1 3) 


0.0114 


0.0114 


0.0003 


(1-13) 


0.0114 


0.0114 


0.0003 


(-1 -1 3) 


0.0114 


-0.0114 


0.0003 


(113) 


0.0114 


-0.0114 


0.0003 


( 4) 


0.0121 


-0.0000 


0.0001 


( 3 -3 -1) 


-0.0075 


-0.0075 


0.0001 


(-3 3 -1) 


-0.0075 


-0.0075 


0.0001 


(3 3-1) 


-0.0075 


0.0075 


0.0001 


(-3 -3 -1) 


-0.0075 


0.0075 


0.0001 



TABLE I: Plane wave expansion coefficients c^(G) for jj, — 
+z. The integers in the first column give Gi, Gr2, Cr3, re- 
spectively the (x, y, z) G cartesian coordinates in units of 
(27r/a). The real and imaginary parts of c+-(G) are displayed 
in columns 2 and 3 respectively. Column 4 shows |c+2(G)p, 
and only coefficients \c+z \ > 10~^ are shown. The coefficients 
for all other minima may be obtained using the symmetry re- 
lations c,i(G) = c!l^(-G); Cx(Gi,G2,G3) = Cz{G3,G2,Gi), 
and Cy{Gi,G2, G3) = Cz{Gi, G3, G2). 



changing the species of the atomic constituents abruptly, 
i.e. across one monolayer (IML). However, the envelope 
function equation allows a continuous choice of the inter- 
face position z/, even within 1 ML width, as discussed 
in Ref. [l^. Although zj is ill-defined within a ML length 
scale, it is convenient to keep this simple interface model 
for it allows decomposing the integral in equation ^ into 
terms that arc easily identifiable, some of which arc famil- 
iar from previous theoretical treatments. Other models 
of the interface have been studied in Ref. [Tol . 

In equation ^ we integrate by parts the term propor- 
tional to U{z) = Ustep{z) given in ([7]), 



- r r pOO 

+ ^ re(.-..)^Me'^M. 

Q J-00 dz 

\^^[z)\^ e'^- zdz 
-00 ^{z) 

= ^U,\^{zj)\'e^^^' (8) 

^^^o«^e-d. (9) 
J Q dz 

\m{z)\^ e'Q'—zAz , (10) 
-00 e(^) 



where we arc interested in an interface bound state, so 
that the electronic density |5'(z)|^ vanishes aX z ^ ±00. 
Three terms are left, labeled ([1]), ^ and (HO]) above. 
Term ([5]) evidences the role of the electronic density at 
the interface |5'(z/)| through a S function at z = z/; 
term © gives the contribution of the evanescent tail of 
the electronic envelope function into the barrier material 
z > z/; term (jlOp represents an intervalley scattering 
induced directly by the electric field, which we find to 
be vanishingly small. Summation of these contributions 
over the reciprocal lattice vectors [see equation ([5])] leads 
to three contributions to the intervalley coupling 



Vvo^Vs + Ve + Vf, 



(11) 



We take for U{z) the step potential, which is the most 
favorable model for the conduction band profile between 
Si and the barrier in terms of maximizing the valley cou- 
pling.^ It is written as^ 

C/(z) = t/,tcp(2) = C/oe(z-z/), (7) 

where z/ is the position of the interface and Uq is the 
conduction band offset. The step potential aims at mod- 
eling a perfectly sharp interface, a concept that involves 



that is, the delta- function contribution Vs, the evanescent 
term contribution Ve, and the electric field contribution 
Vf. 

Besides the mismatch between the conduction band 
minima of the two materials, we also consider the change 
in dielectric screening constant from the semiconductor 
to the insulator. The dielectric function is taken as 
e(z < z/) = esi and e(z > z/) = ebarrior- Note that 
e(z) also introduces a kink in the electrostatic potential 
in the scale of the monolayer separation, which could in 
principle contribute to the intervalley coupling. 



III. NUMERICAL SOLUTION AND RESULTS 



Envelope Function Variational Approaches: 
Finite Differences and Trial Function 



The envelope function "^{z) is obtained from equa- 
tion which has no analytic solution. The approach 
adopted to solve this equation must be carefully chosen, 
since the intervalley coupling depends explicitly on the 
envelope function details, as shown in equations (|8]), (|9]) 
and (fTO]). The first term (jS]), which gives rise to Vs in 
equation (|lip . is proportional to the electronic proba- 
bility density at the interface , while the second 
term ^ depends on the envelope tail inside the barrier, 
defining the contribution Ve of the evanescent tail. The 
term ()10|) is always found to be negligibly small compared 
to the others, and is not included in the results presented 
here. 

In order to get an analytic approximation for the en- 
velope function, we present initially results obtained us- 
ing the variational approach; we tried several functional 
forms for the variational envelope function, all satisfying 
the following boundary conditions: in the far semicon- 
ductor region (z — > — oo), under a constant electric field, 
the envelope function is approximated by a gaussian de- 
cay; towards the barrier material the short range decay 
{z > z/) is nearly exponential. The trial function that 
typically gave the lowest energy expectation value was 



,z < Zl 

,Z > Zl 



(12) 

with N, = {Jdz *|) (i=A,B). 

If we take the barrier to be infinite, the optimized vari- 
ational parameters become zq = zj and B = 0. But in 
general the barriers are finite, some penetration is ex- 
pected so that Zq > zj and B > 0. From the continuity 
conditions '^a{zi) = ^'^(z/) and ^'^(z/) = ^'^(z/), and 
the normalization + B^ = 1, we obtain expressions 
for A, B, and zq in terms of the variational parameters 
a and /3; e.g. For zi = 0, we have 



zo 



4a 



(13) 



A = 




8z2 a3/2 



P\ 4zoVa-Fe2^o"v^erfc(V2azo) 



(14) 



B 



Izl a3/2 



4zoVa(2zoa + P) + e2^o"^/3erfc(V2azo) 
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FIG. 1: (Color online) Energies and envelope functions for 
a conduction band offset C/q = 3 eV. (a) Comparison be- 
tween the expectation value for the ground state energies ob- 
tained variationally by solving the Finite Differences equation 
through the steepest descent method (squares) and using the 
trial wavefunction in equation H12p (diamonds). Exact results 
for (7o — > oo from Ref . [stI. are also shown (solid line), (b)-(d) 
Envelope functions squared (probability densities) for the in- 
dicated electric fields. The solid line is obtained solving the 
Finite Differences equation variationally (through the steep- 
est descent method). The dashed line is obtained from the 
trial-function in equation (|12p . In (b), F — 1.00 V/nm and 
the two approaches are in good agreement, while at lower elec- 
tric fields, (c) and (d), the functional form in equation (|12p 
leads to somewhat more localized states near the interface. 



where erfc is the complementary error function. 



2 

rfc(a;) = 1 - erf(a;) = 1 - — / . 



dt. 



(16) 



(15) 



Analytic solutions for the case of z/ 7^ may also be 
obtained. However, for the step function potential and 
at relatively large length scales, the interface position 
is irrelevant: it separates two semi-infinite regions, one 
filled with Si and the other with the barrier material. At 
atomic distance length scales, the coupling dependence 
on z/ shows oscillatory behavior with a period of 1 ML, 
a peculiarity of the EMA combined with the underlying 
Bloch states, as discussed in detail in Ref. [l^. Data pre- 
sented there (Figure 5 of Ref. [l^ indicate that, for the 
step potential, these oscillations cause an uncertainty of 
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about 20% in IVyol- In what follows, for simplicity and 
definiteness, we arbitrarily fix the value of z/ = bearing 
in mind that our results for the coupling are not accurate 
within up to ^ 20%. 

In order to provide a more robust estimate of the wave- 
function, independent of our choice of its functional form, 
we also solve equation ([2]) through a finite differences 
method. This approach allows us to obtain the wave 
function variationally without inputing a guess for the 
functional form of the envelope. The envelope function 
is discretized and the derivatives are approximated by 
the slope of the interpolated straight line between two 
discrete points. Each discrete point is taken as a sepa- 
rate variational parameter, under the constraint that the 
wave function is normalized. The minimum energy in 
this configuration space is obtained through the Steep- 
est Descent method. This strategy should in principle 
permit us to discuss a wide range of electric fields. 

The trial function and the finite differences expecta- 
tion value of the energy are given in Fig. [TJa), according 
to which the lowest value down to F ^ 0.01 V/nm is ob- 
tained within finite differences/steepest descent method. 
The exact result for infinite barrier, given in Ref. [stI . is 
also included, and a good agreement with our numeri- 
cal results is obtained. As the field increases and pushes 
the wavefunction towards the barrier, an evanescent tail 
into the barrier region is formed for C/g = 3 eV, lowering 
the electronic energy with respect to the impenetrable 
barrier C/q — > oo case, as expected [upper right data in 
Fig-HHa)]. The good agreement for the expectation value 
of the energy among the three methods, for example at 
F=0.1 V/nm in Fig 1(a), does not imply that the wave- 
functions are in agreement - as illustrated in Fig [Ijc) 
and as is obvious for the infinite barrier case, where the 
wavefunction for z > is exactly zero. It is known that 
the energy alone is not a valid criterion to guarantee the 
wavefunction validity, in particular here for z > - from 
which the coupling is calculated. Imposing a predeter- 
mined form to the envelope most certainly will give spu- 
rious results for Vs and Ve- This points the numerical 
finite differences as the most adequate approach, with 
the additional capability of describing different interface 
profiles 

The finite differences approach is expected to be more 
accurate than trial-function based methods. However, 
in practice, numerical constraints limit its implementa- 
tion and reliability. Below F « 0.01 V/nm (not shown), 
we find that the lowest energy actually corresponds to 
Eq. (jl2p . a consequence of our numerical accuracy lim- 
itations. Comparison between the envelope functions 
squared (probability densities) for a range of electric 
fields and for a conduction band offset Uq = 3 eV is given 
in Fig. mb)-(d). It is clear from the figures that the en- 
velope defined in Eq. is in very good agreement 
with the one obtained from Finite Differences for rela- 
tively high electric fields [this is illustrated in Fig. HKb)], 
but the difference at lower electric fields is noticeable 
[see Fig. [IJc), (d)]. Results presented below correspond 



to the Finite Differences solution to equation (2) in the 
F > 0.01 V/nm regime. 



B. Contribution from Vf, Vs and Ve 

We now proceed to calculate the three terms in equa- 
tion (|lip . As mentioned above, we find the electric field 
contribution Vp to be orders of magnitude smaller than 
the other two terms for any value of ebarrier- This is 
an indication that the kink in the electrostatic potential 
introduced by the change in dielectric constants is not 
singular enough to produce sizable coupling between the 
±z valleys. 

The other two terms arc also left unchanged if the same 
conduction band offset Uq is imposed and different val- 
ues of ebarrier are adopted, meaning that the electronic 
wavefunction in all cases does not penetrate the barrier 
material deep enough to be affected by the electric field 
inside it. Therefore, we disregard the electric field in the 
barrier. This is beneficial to our model, since we can 
characterize the barrier material by the conduction band 
offset alone, and not involve other properties specific to 
the barrier material. 

A richer behavior is obtained from V5, the contribution 
from the 5-function given in Eq. ([8]). Since this term is 
proportional to the product Uq \'^{zj)\ , there is a trade- 
off between the barrier height and the envelope function 
penetration into the barrier material. In Fig. [5] we note 
that jV^I increases with Uq, meaning that, for F = 10^ 
V/nm, the increase in Uq prevails over the reduction in 
|*(~^/)|^ 

Finally, the term Ve arising from the evanescent tail 
contribution 5'(z > z/) also presents a non-trivial trade- 
off. While some penetration of the envelope function 
into the barrier material is needed for this term to be 
non- vanishing, the integrand in equation ^ is highly os- 
cillatory so that if the envelope function penetrates more 
than a few monolayers, Ve integrates to zero. Since SiGe 
barriers present fairly large electronic penetration, the 
calculated Ve for these materials is much smaller than 
in the case of Si02 barriers. Besides, in the SiGe barrier 
the Ve contribution only gives a change in the complex 
phase of Vyo, thus leading to \Vvo\ ~ I ^5 1- Still, since 
this small Ve contribution to the intervalley coupling in 
Si/SiGe heterostructurcs changes its complex phase, it 
leads to a different ground state combination of the ±z 
valleys. 



C. Umklapp Processes 

Since we take the plane wave expansion of both ±z 
conduction band minima, Umklapp processes are fully in- 
cluded here. An estimate of their relevance is given by the 
total summation over them, which gives IVajumkiapp = 
0.02907C/o|*(z/)|^ao, while the contribution from G = 
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G' = is IV^Ir = 0.01108 f/o qq. 

We also look at each plane wave contribution sepa- 
rately. The most prominent terms in the PW expansion 
of the Bloch functions at the conduction band minima are 
not the G = term, but the first nearest neighbors in the 
BCC reciprocal lattice, as can be seen in Table HI 

Table III shows the numerical prefactors, multiplying 
L^o iQ, of the contributions to the imaginary part, 

confirming that the most relevant contributions come 
from the F point and its first eight neighbors in the BCC 
reciprocal lattice of Si since we get over 70% of the to- 
tal value coming from these points. Notice the Umklapp 
terms are the most relevant. 

It is also interesting to note that for the imaginary 
part the most relevant contributions are obtained for 
G = G'. This is equivalent to taking the zeroth or- 
der plane wave expansion of the product it!j_(r)u_(r), a 
commonly adopted approximation in the EMA theory of 
shallow donorsi^ But a complete analysis of the Umk- 
lapp processes reveals that this is particular to the delta 
contribution, and that there is no general justification 
for this approximation for the interface induced valley 
coupling. 



G(27r/a) 


G'(27r/a) 


Contribution ((7o ^(^-r)^ o-o) 


(0,0,0) 


(0,0,0) 


-0.01108 


(1,1,1) 


(1,1,1) 


0.00410 


(-1,1,1) 


(-1,1,1) 


0.00410 


(1,-1,1) 


(1,-1,1) 


0.00410 


(1,1,-1) 


(1,1,-1) 


0.00410 


(-1,-1,1) 


(-1,-1,1) 


0.00410 


(-1,1,-1) 


(-1,1,-1) 


0.00410 


(1,-1,-1) 


(1,-1,-1) 


0.00410 


(-1,-1,-1) 


(-1,-1,-1) 


0.00410 



TABLE II: Most relevant contributions to the imaginary part 
of Vg. The sum of these terms is 0.02175 Uo \'^(zi)f ao, un- 
derestimating this term by only 25%. The real part of this 
contribution vanishes due to our choice of coordinate system. 



D. Complete Valley Coupling 

The absolute values, \Vs\ and IVeI, of the two terms 
contributing to VVo and \Vvo\ — \Vs + Ve\ arc shown 
as a function of the barrier height in Fig. [51 It is clear 
here that in general the two terms have different behav- 
iors, and that they both give important contributions to 
the total coupling. However, for estimating the valley 
splitting, |Vvo| = I ^5 1 is a reasonable approximation for 
small barriers, such as those in Si/SiGe heterostructures. 
Noting that \Vvo\ < iKs] -I- \ Ve\, we infer that Vyo is in 
general a complex number, not a purely imaginary quan- 
titj^. Since Vs is purely imaginary, this triangle inequal- 
ity is only true if Ve has a non- vanishing real part. For 
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FIG. 2: (Color online) Intervalley coupling as a function of the 
conduction band offset Uo- All data correspond to an external 
field of 10^^ V/nm. The absolute value of the delta function 
contribution \ Vs \ and the evanescent tail contribution \ Ve\ are 
depicted separately, as well as the real and imaginary parts 
of the total coupling Wo = Vs + Ve- The offset for Si/SiGe 
interfaces is 0.1 — 0.2 eV and for Si/Si02 interfaces is 3 eV. 



instance, at Uq = 3eV, we have an intervalley coupling of 
Wo = (-0.100 + i0.158)meV. Also, Vs and Ve increase 
monotonically with Uq , while Wo decreases at large off- 
sets. This indicates that the relative phase between Vg 
and Ve changes with Uq- 

In principle, SiCe and Si02 barriers could lead to sim- 
ilar intervalley couplings. Of course these two materi- 
als present very different interface morphologies and are 
grown with different techniques, which should lead to dif- 
ferences between the intervalley couplings measured in 
each design. There has been reports of significantly dif- 
ferent valley splitting even for the two interfaces of the 
same quantum welli^ 

The valley coupling calculated for a range of external 
electric fields, covering several orders of magnitude are 
summarized in Fig. [3] This range expands by two orders 
of magnitude the previously reported range in Ref. [l^ 
(see inset). We take the upper limit of F just below the 
Si02 breakdown field F ^ 3 V/nm. This bound is indi- 
cated in the figure, as well as the SiGe breakdown field 
w 10~^ V/nm. The lower value of F is constrained by 
our numerical accuracy, as discussed in Sec. IIII Al For 
large enough values of F, the valley splitting can be of 
the order of 10 meV, which is compatible (within error) 
with the giant splitting observed by Takashina et al^ 
This means that the results obtained in a SIMOX inter- 
face could be related to nearly perfectly sharp interfaces 
combined with relatively high gate fields, which might be 
consistent with the experimental conditions in terms of 
these two parameters. Other effects may be present that 
can further enhance this coupling and account for or con- 
tribute to very large splittings, as the recently proposed 
mechanism involving interface states <^ 

Previous studie a^^'^^ also report a linear dependence 
\Vvo\ = AF, where A is a model-dependent length. This 
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behavior in difFercnt models is qualitatively and quanti- 
tatively addressed in Sec. IIVI 
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FIG. 3: (Color online) Data points give the calculated inter- 
valley coupling as a function of applied electric field, both 
in a log scale, for barrier height Uq = 3 eV. A hnear fit for 
the data points with Wo given in meV and F in V/nm, is 
included (dotted line through the points), leading to an es- 
timated A ~ O.27A. The vertical dashed lines correspond to 
the experimental values of the breakdown field for SiGe and 
Si02. The dotted horizontal line represents the measured 
valley coupling (or half the splitting) reported in Refl2ll. The 
solid straight line is a possible prediction from the model in 
Ref. [2^, calculated from equation (|17p. taking a — 0.43 A 
and = F- The inset shows the same results presented 

in Ref. [13, covering a much narrower range of field values, for 
IJa = ?> eV (solid line) and 125 meV (dotted line) 



IV. DISCUSSIONS AND CONNECTIONS WITH 
PREVIOUS MODEL CALCULATIONS 

The problem of interface-induced intervalley coupling 
has been treated extensively in the literature, particu- 
larly within EM A. Since the complete problem is not ex- 
actly solvable, several approximations regarding the bar- 
rier potential and the nature of the electronic states have 
been adopted in different contexts. We discuss in this 
Section connections between our study and some previ- 
ous contributions. 



A. Electric Field Dependence of the Valley 
Splitting 

The applied field F is one of the key parameters af- 
fecting VvOi and probably the most controllable. Early 
work by Sham and Nakayama^^ established a linear de- 
pendence of the valley splittting on F . They address 
the Si valley splitting problem from the perspective of 
the electrons in a space-charge layer in a MOSFET. The 



effect of the Si-Si02 interface on electron dynamics is 
studied by considering incoming (towards the interface) 
and outgoing (away from the interface) Bloch waves in 
the Si region and evanescent waves decaying from the 
barrier inside the semiconductor. The potential barrier 
is modeled to be infinite, which disregards the evanes- 
cent tail of the electronic wavefunction into the barrier. 
Also, their study is devoted to the case of a planar den- 
sity of electrons at the interface (of the order of 10^^ 
cm~^ electrons), including many body corrections. In 
contrast, we develop here a theory for a single electron in 
a bound state, relevant in the context of quantum com- 
puting. Therefore, any comparison with Ref. [2^ should 
be taken with caution. 

The valley splitting obtained by Sham and 
Nakayama^^ 



AE = 2\Vvo\ 



(17) 



is proportional to the mean value of the derivative of the 
self-consistent potential. The parameter a, a character- 
istic length related to the intervalley scattering matrix, is 
given as a function of the interface position with respect 
to a crystal (001) plane, (called zq there), showing sharp 
variations as zg runs over a 1 ML range. But the concept 
of an interface as an abrupt change in conduction band 
energy is not well defined in this scale. This ambigu- 
ity comes from the continuum EMA approach combined 
with the atomistic band structure input. These same 
ingredients affect our results, with the calculated valley 
coupling oscillating as a function of z/J^ We estimate 
that the mean value of the derivative of the self-consistent 
potential should be of the order of the field, (^) = F, 
and take \a\ — 0.43 A as estimated in Ref. [2^ Therefore 
\{SkN) = a/2 = 0.21 A. 

More recent EMA-based models^^i^^i^"^ rely on an ef- 
fective coupling potential responsible for the intervalley 
coupling, where the perturbation potential is taken as a 
(5-function with strength obtained from either TB atom- 
istic calculations or experimental data. In our formalism, 
it means taking in equation (|lip only the term V5. As 
discussed in the previous section and shown in Fig. [21 
this approximation is better for the relatively low SiGe 
barriers, but for barrier materials with higher conduc- 
tion band offsets, sizeable contributions (actually reduc- 
ing I VVql) come from the Ve term. The valley couplings 
in Refs. [6l| andjo^ are also found to be linear with elec- 
tric field, and their estimated values for A are 0.72 A and 
1.36 A respectively. 

In our model calculations, the direct {kz,k^z\ cou- 
pling {Vp) mediated by the field term in H is found to 
be negligible. The effectiveness of F comes from the elec- 
tronic charge being pulled towards the barrier material, 
increasing the wavefunction amplitude at and beyond the 
interface position. These shifts contribute to the terms Vs 
and Ve in ((TT|) . From the linear fit of our data in Fig. |31 
we get A (present wor k)=27.4 mV/(V/nm) = 0.27A, thus 
comparable to \{Sk^N)^ Results in Ref. |6l| (alterna- 
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tively, Ref. [63 ) are larger than SkN by a factor of 3 (6) , 
and are 2 (5) times larger than ours. Of course, there 
are many differences between the models which may ac- 
count for the different results, but it is interesting to note 
from Fig. [2] that the delta- function term alone, \Vs\, al- 
ways overestimates the full coupling \Vvo\ = \Vs + Ve\- 
Therefore, considering only the delta function potential 
in the coupling term may be related to larger values pre- 
dicted for A in Refs. HH and HI. 

The TB approach is also used to obtain more accu- 
rate boundary conditions at the interface for the EMA 
wavefunctioui^ Furthermore, the ambiguity introduced 
by the interface position discussed before was also rec- 
ognized in this study, and explored in Ref. HH to treat 
the interface position within the uncertainty range as a 
fitting parameter to match the TB and effective mass 
wavefunctions in a finite quantum well. 



B. Atomistic Approaches 

More realistic quantitative estimates may be obtained 
from atomistic many band TB description of Si and the 
barrier material;^"— which may also account for inter- 
face disorder such as alloying effects and surface rough- 
ness. Such approaches are not always adequate or meant 
to give a general picture of the physical mechanisms be- 
hind the intervalley coupling, since the focus is to achieve 
quantitative accuracy. In some cases reasonable trade-off 
between atomistic description and analytical interpreta- 
tions has emerged from simple one dimensional two band 
TB models^iia 

Usually, the TB simulations are performed within a su- 
percell approach, which involves periodic boundary con- 
ditions. So a Si/barrier interface is actually modeled by 
a finite-width Si slab surrounded laterally by two barri- 
ers, implying a quantum well arrangement if the barri- 
ers are wide enough (or superlattice for thin barriers). 
This introduces a width parameter to the Si slab which 
was shown to play an important role in the intervalley 
coupling except for high external electric fields and 
wide enough Si slabs, so that the electron interacts with 
a single interface. On the other hand, within EMA it is 
possible to pinpoint the role of a single interface, as both 
Si and the barrier material correspond geometrically to 
semi-infinite slabs. 

Despite these differences, some TB studies give em- 
pirical evidence about the EMA analytical insights ob- 
tained here: (i) the linear behavior of the valley splitting 
with electric field, as discussed in Sec. IIV Al (ii) the rela- 
tionship between the penetration of the wave-function in 
the barrier material and the barrier splitting explained in 
Sec. lIIII and quantified by Eqs dS]) and © was also empir- 
ically obtained by Srinivasan et al^ for Si electrostatic 
quantum dots embedded in SiGe buffers; (iii) finally the 
connection between the conduction band offset and the 
valley splitting demonstrated here in Fig. [5] similar to the 
TB study reported in Ref. [3^ 



These examples illustrate how EMA studies may con- 
tribute to clarify the physics behind phenomena numer- 
ically obtained by more detailed atomistic approaches. 

Other atomistic studies include ingredients that were 
not implemented in our formalism, such as magnetic field 
and interface disorder, so that it is not possible to quan- 
titatively compare the results demonstrated within our 
general EMA study to those obtained within TB for par- 
ticular systems/geometries. 

C. Low Field Dependence of the Valley Splitting 

A non-linear dependence of VVo with F has been pre- 
sented in TB studies by Grosso et al^ and Boykin et 
al^ at low fields. We can not assess this behavior be- 
cause of instabilities in our numerical procedure, mainly 
due to the size of the simulation cell required by the wave- 
function spread at very low electric fields. It is possible 
that this nonlinear behavior is connected to the quan- 
tum well behavior at low F, since the models in Refs. [33l 
and refer to a quantum well, involving lengths as- 
sociated with the Si well and the barriers widths, on 
which (particularly the Si well width) the TB results are 
quite sensitive. In all cases a linear behavior is obtained 
at l arg er F values. The linear behavior is predicted in 
Refiifrom Eq. obtained there, with (|f) = F. It 
also emerges from a model for the variational envelope 
in a triangular-type potential with infinite barrier poten- 
tial, simpler but still similar to Eq. (fT2|). proposed by 
Fang-Howard^ and generalized by Friesen et al^ to al- 
low some penetration probability into the barrier region. 
In this case an analytic solution is obtained, and it can be 
shown that |5'(2/)p is proportional to the electric field, 
leading to linear increase of the intervalley coupling with 
F. However, in most cases of interest the applied electric 
field should lead to a large enough VVo 7 thus well within 
the linear regime {F > 10"'^ V/nm here). 

V. SUMMARY AND CONCLUSIONS 

We presented an EMA-based study of the valley split- 
ting in Si induced by a Si/barrier interface. Our ap- 
proach combines EMA with calculated Bloch functions, 
for which we give values of the relevant plane-wave ex- 
pansion coefficients. The range of splittings 2 1 VVo I we 
obtained are in fair agreement with measurements in 
Si/Si02 and Si/SiGc interfaces32i Our results for \Vyo\ 
arc comparable with experimentally measured values, in- 
dicating that we have probably included the most rel- 
evant physical ingredients in our model, as discussed 
in Ref. and expanded here. In particular, we show 
here that the puzzling values of the valley splittings ob- 
tained in Ref. |2ll could result from particularly sharp 
interfaces in the very high electric field regime. We also 
confirm the linear dependence of the valley splitting on 
the electric field in this regime. Nonetheless, many ef- 
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fects contributing to the valley splitting are not included, 
such as strain, interface misorientation;^ atomic scale 
disorder,— lateral confinement,— or many-body correc- 
tion o^^'^^ and the recently proposed contribution from 
interface states 

It is important to reiterate here the double-focus of 
this paper : (i) to examine the effect of an applied elec- 
tric field on the valley-orbit coupling at the interface, 
and (ii) to provide a clear interpretation of the physics of 
valley splitting. Understanding these points would hope- 
fully reveal fundamental elements that affect the electron 
valley coupling and guide the identification of a suitable 
environment for electron spin qubits. By using the ef- 
fective mass approximation we implicitly assume that 
the atomic configuration at the interface is not changed 
by the applied field (which should be a good approx- 
imation at low fields), and that interface roughness is 
small at the length scale over which the electron wave 
function changes significantlyi^ A comprehensive study of 
valley splitting, including the effects of interface rough- 
ness, thickness, applied electric field, and other elements 
mentioned in the previous paragraph, would necessarily 
contain an atomistic component. Such a study would 
be required to clearly establish the limit of applicability 
of EMA, which is essentially a mean field approach for 
the bound electron wave function, interface disorder and 
composition profiles. 

Our EMA approach provides a clear interpretation of 
the physics of valley splitting, revealing fundamental el- 
ements affecting the coupling and eventually providing a 
suitable environment for electronic spin qubit operation. 
While the investigation of other effects mentioned above 



is desirable, a simpler model that goes beyond the phe- 
nomenological approach, as provided here, is useful in 
guiding nanofabrication and device operation efforts. In 
particular, a profound understanding of the valley cou- 
pling induced by interfaces may pave the road to the 
quantum manipulation and processing of the valley de- 
gree of frccdomi^ 

We also provide a critical analysis of the valley physics 
theory available in the literature, and discuss some of 
the ingredients that are imperative in a comprehensive 
theory. These include the correct interpretation of the 
valley coupling as a complex number and the inclusion 
of Umklapp processes. 

In conclusion, we have calculated electron valley split- 
ting in Si at a Si/barrier interface. We show that a size- 
able single-particle valley-orbit coupling can be obtained 
by applying a high enough external field and choosing an 
optimal barrier material providing a suitable potential 
barrier height and high quality abrupt interface. 
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